Given a Gaussian oscillator of the form:
$$\epsilon_2 = A\left[e^{-4\ln(2)\left(\frac{\omega-\omega_0}{\sigma}\right)^2} - e^{-4\ln(2)\left(\frac{\omega+\omega_0}{\sigma}\right)^2}\right]$$as defined in [1] so that $\epsilon(\omega)$ is Kramers-Kronig consistent by calculating $\epsilon_1$ with a KK analysis, we want to obtain the spectral weight of such an oscillator. The spectral weight is given by
$$SW_{G} = \int\limits_0^{\infty}\sigma_1(\omega) d\omega = \int\limits_0^{\infty}\omega\epsilon_0\epsilon_2(\omega) d\omega $$where $\sigma_1$ is the real part of the optical conductivity and $\epsilon_2$ is the imaginary part of the dielectric function.
from sympy import *
init_printing()
omega = symbols('omega')
A, sigma, omega_0 = symbols('A sigma omega_0', real = True, positive = True)
First, to get familiar with the integrals, the strategy will be to calculate $\int\limits_0^{\infty}Ae^{-C\left(\frac{\omega\mp\omega_0}{\sigma}\right)^2}$ with $C = 4\ln(2)$.
C = symbols('C', real = True, positive = True)
I_m = Integral(A*exp(-C*((omega-omega_0)/sigma)**2), (omega, 0, oo))
I_p = Integral(A*exp(-C*((omega+omega_0)/sigma)**2), (omega, 0, oo))
I_m-I_p
Performing a change of variables $x = \frac{\omega\mp\omega_0}{\sigma}$
x = symbols('x')
I_m = I_m.transform((omega-omega_0)/sigma, x)
I_p = I_p.transform((omega+omega_0)/sigma, x)
I_m-I_p
The area of integration gets reduced to the interval $[-\frac{\omega_0}{\sigma}, \frac{\omega_0}{\sigma}]$. And exploiting the even integrand, it ends like
2*A*sigma*Integral(exp(-C*x**2), (x, 0, omega_0/sigma))
Performing the integral
_.doit()
Now, performing the real integral for the Gaussian spectral weight. Again splitting the integral in two pieces with $C = 4\ln(2)$.
epsilon_0 = symbols('epsilon_0', real = True, positive = True)
NumCoeff = A*epsilon_0*sigma
SW_m = Integral(NumCoeff*exp(-C*((omega-omega_0)/sigma)**2), (omega, 0, oo))
SW_p = Integral(NumCoeff*exp(-C*((omega+omega_0)/sigma)**2), (omega, 0, oo))
SW_m-SW_p
Performing the same change of variables as before $x = \frac{\omega\mp\omega_0}{\sigma}$
SW_m = SW_m.transform((omega-omega_0)/sigma, x)
SW_p = SW_p.transform((omega+omega_0)/sigma, x)
SW_m-SW_p
For the term linear in $x$ (odd), the interval of integration is again $[-\frac{\omega_0}{\sigma}, \frac{\omega_0}{\sigma}]$, which renders the integral to $0$. Reordering the rest of the integrals and keeping the leading $A\epsilon_0\sigma$ coefficient out.
SW_1 = Integral(omega_0*exp(-C*x**2), (x, -omega_0/sigma, oo))
SW_2 = Integral(omega_0*exp(-C*x**2), (x, omega_0/sigma, oo))
SW = SW_1 + SW_2
SW
This is the same integral with overlapping intervals. Once in the range $[-\frac{\omega_0}{\sigma}, \frac{\omega_0}{\sigma}]$ and twice in the range $[\frac{\omega_0}{\sigma}, \infty]$. Absorbing $\omega_0$ to the numerical prefactor and leaving it out for clarity.
NumCoeff *= omega_0
SW_1 = Integral(exp(-C*x**2), (x, -omega_0/sigma, omega_0/sigma))
SW_2 = Integral(exp(-C*x**2), (x, omega_0/sigma, oo))
SW = SW_1 + 2*SW_2
SW
The first integral is already known from the preparatory calculations, but we recalculate it again here. This give us a closed form for the spectral weight of a Gaussian oscillator.
SW = SW.subs(C, 4*ln(2))
result = SW.doit().simplify()
result *= NumCoeff
result
[1] De Sousa Meneses, D., Malki, M., & Echegut, P. (2006). Structure and lattice dynamics of binary lead silicate glasses investigated by infrared spectroscopy. Journal of Non-Crystalline Solids, 352(8), 769–776. doi:10.1016/j.jnoncrysol.2006.02.004
Produced by Ignacio Vergara Kausel vergara at ph2.uni-koeln.de
Available at https://github.com/ivergara/science_notebooks
Licensed under GPL3